function [lambda, gamma]=rbf(S,Y,flag)
% RBF插值模型
% INPUTS:
%   S    - 样本点矩阵 (m x d)，其中 m=样本点数, d=维度
%   Y    - 样本点对应的目标函数值 (m x 1)
%   flag - 指定使用的 RBF 类型 (字符串)，可选值：
%          'cubic'        - 三次径向基函数
%          'TPS'          - 薄板样条 (Thin Plate Spline)
%          'linear'       - 线性径向基函数
%          'Gaussian'     - 高斯径向基函数
%          'multiquad'    - 多二次径向基函数 (Multiquadric)
%          'invmultiquad' - 反多二次径向基函数 (Inverse Multiquadric)
%
% OUTPUTS:
%   lambda - 径向基函数权重向量 (m x 1)
%   gamma  - 多项式项系数向量 ((d+1) x 1)
%
%--------------------------------------------------------------------------

[m,n]=size(S);
P=[S,ones(m,1)];

R=zeros(m,m);
% 使用欧几里得距离矩阵 (优先调用 pdist2)
if exist('pdist2', 'file')
    R = pdist2(S,S,'euclidean');
else
    % 如果没有 pdist2，则手动计算距离
    for i = 1:m
        for j = 1:m
            R(i,j) = sqrt(sum((S(i,:) - S(j,:)).^2));
        end
    end
end

% 根据不同的 flag 选择 RBF 核函数
if strcmp(flag,'cubic') % 三次径向基
    Phi= R.^3;
elseif strcmp(flag,'TPS') % 薄板样条
    R(R==0)=1; % 避免 log(0) 出现 NaN
    Phi=R.^2.*log(R);
elseif strcmp(flag,'linear') % 线性径向基
    Phi=R;
elseif strcmp(flag,'Gaussian') % 高斯径向基
    rho = max(max(R))/(n*m)^(1/n);
    Phi=exp(-R.^2/rho^2);
elseif strcmp(flag,'multiquad') % 多二次径向基
    rho = max(max(R))/(n*m)^(1/n);
    Phi=sqrt((R.^2 + rho^2));
elseif strcmp(flag,'invmultiquad') % 反多二次径向基
    rho = max(max(R))/(n*m)^(1/n);
    Phi=1./sqrt(R.^2 + rho^2);
end

% 构造线性方程组 A * params = RHS
A=[Phi,P;P', zeros(n+1,n+1)];
RHS=[Y;zeros(n+1,1)];

% 数值稳定性改进：检查矩阵条件数，必要时加正则化
cond_num = cond(A);
if cond_num > 1e12
    reg_param = 1e-8;
    A_reg = A + reg_param * eye(size(A));
    params = A_reg\RHS;
else
    params = A\RHS;
end

lambda=params(1:m);
gamma=params(m+1:end);

end % function
